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The solidification of polycrystalline materials can be modelled by orientation-field models, which 
are formulated in terms of two continuous fields: a phase field that describes the thermodynamic 
state and an orientation field that indicates the local direction of the crystallographic axes. The 
free-energy functionals of existing models generally contain a term proportional to the modulus 
of the orientation gradient, which complicates their mathematical analysis and induces artificial 
long-range interactions between grain boundaries. We present an alternative model, in which only 
the square of the orientation gradient appears, but in which the phase and orientation fields are 
coupled by a singular function that diverges in the solid phase. We show that this model exhibits 
stable grain boundaries whose interactions decay exponentially with their distance. Furthermore, 
we demonstrate that the anisotropy of the surface energy can be included while preserving the 
variational structure of the model. Illustrative numerical simulations of two-dimensional examples 
are also presented. 



I. INTRODUCTION 

In polycristalline materials that consist of multiple 
crystalline grains of the same thermodynamic phase, the 
grain size and texture of the material largely determine 
mechanical properties such as toughness and yield stress. 
Therefore, a large effort has been devoted to understand- 
ing the relation between the processing conditions during 
solidification and subsequent ageing of a material and the 
resulting grain structure. From a fundamental viewpoint, 
this is an example of self-organization [l[ : a complex 
structure emerges from a structureless initial state (the 
liquid), and the resulting network of grain boundaries 
evolves with time and coarsens 0, Q . 

Despite a large body of work, many questions regard- 
ing this self-organization process remain open. As often 
in such complex problems, a comparison between exper- 
iments and numerical modeling could be highly useful. 
Therefore, it is desirable to develop quantitative mod- 
eling tools for the solidification and evolution of poly- 
cristalline materials. Among the many modelling ap- 
proaches, the phase-field method has rapidly gained in 
popularity in recent years due to its versatility and ro- 
bustness. For recent reviews on its applications in mate- 
rials science, see Refs. 0-0] • 

Numerous phase-field models for the description of 
polycristalline solidification and grain growth have al- 
ready been developed. They can be grouped in two cate- 
gories. In the so-called multi-phase-field approach 0411!], 
each grain is described by an individual phase field. The 
grain boundary energies are determined by the coupling 
between the phase fields. Despite the necessity to handle 
a large number of phase fields, this approach can be made 
quite efficient from a numerical point of view [ToL [T^ . 
Nevertheless, it suffers from the fact that for each simu- 
lation the coupling coefficients between the phase fields 
need to be computed (and stored) in order to reproduce 
the proper grain boundary energies Moreover, since 
the crystalline orientation in each grain is fixed, phenom- 



ena such as grain rotation and plastic deformation cannot 
easily be described. Finally, from a fundamental point of 
view, a large number of fields certainly does not consti- 
tute the most "economic" description of a polycristalline 
structure, since the state of matter in any space point 
can be specified by a set of a few variables, such as the 
degree of order and the local lattice orientation. 

As a consequence of the latter consideration, alterna- 
tive models have been developed that use a single phase 
field coupled to an orientation field p^4f8| . The phase 
field describes the local thermodynamic state - either 
solid (<fi — 1) or liquid (0 = 0)- and the orientation field 
is used to describe the local orientation of the crystallo- 
graphic axes in the solid with respect to some fixed ref- 
erence system. In two dimensions, the orientation field 
is simply an angle 8{x.,i), while in three dimension a 
more complex formalism is needed [l9l [20j . In this frame- 
work, a grain boundary is a diffuse but localized region in 
space between two regions of the same thermodynamic 
state (both are solid) but different crystalline orienta- 
tions. The change of the crystalline orientation is con- 
centrated in the grain boundary where the phase field is 
no longer equal to 1. 

While this kind of description is appealing, it turns out 
that the construction of a viable model is far from simple 
because of the need to obtain localized grain boundaries. 
Let us consider the two-dimensional (2D) case with a 
scalar angle field and try to construct a free energy func- 
tional that leads to localized solid-liquid interfaces and 
grain boundaries. For the solid-liquid interfaces, this is 
straightforward: since there are two privileged and dis- 
tinct thermodynamic states (solid and liquid) , diffuse in- 
terfaces arise from the free energy balance between a 
double- well potential and a square-gradient term for the 
phase field. In the case of grain boundaries, matters are 
not so simple. Indeed, there is no privileged angle and 
the free energy must be rotationally invariant. Therefore, 
a potential that depends on the angle is forbidden, and 
only terms that contain spatial derivatives of 9 can be 
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used in the free energy. Then, in order to model local- 
ized grain boundaries, at a first glance it seems sufficient 
to introduce a term that couples the energy cost of an- 
gular variations to the local value of the phase field and 
induces a sufficiently high energy penalty for angle vari- 
ations in the solid. 



However, a simple scaling argument, reviewed in de- 
tail below, shows that localized grain boundaries cannot 
be obtained by a straightforward coupling of a square- 
gradient term in 9 to the phase field. In fact, such a model 
is equivalent to the usual thermodynamic description of 
nematic liquid crystals (2lj . in which indeed no "grain 
boundaries" (localized changes in the director field) exist. 
The solution proposed in the previous orientation-field 
models is to include an additional term proportional to 
|V#| in the free energy. Localized grain boundaries then 
arise from the balance between this term and the square 
gradient in 9. While this approach yields localized grain 
boundaries, it also has several problems. First, from a 
fundamental viewpoint it seems difficult to justify such 
a form in the spirit of a Landau expansion. Second, the 
functional derivative of this non-analytic term leads to a 
singular diffusion equation for 9 and to non-local inter- 
actions between grain boundaries [22j]. The latter effect 
can be removed by a proper regularization [l6j . but this 
introduces additional parameters into the model. 

Here, we propose an alternative approach. Our model 
contains only a square gradient term in 9, but relies on 
the use of a singular coupling function between orienta- 
tion and phase field to make the energy cost of angular 
variations in the solid effectively infinite. Similar ideas 
were proposed some years ago by Lusk [23j, but were to 
our best knowledge not pursued any further until now. 
The interesting point about this approach is that the 
equations for the equilibrium grain boundaries and solid- 
liquid interfaces are ordinary (non-singular) differential 
equations that can be analyzed using standard mathe- 
matical tools and concepts from dynamical systems the- 
ory. Furthermore, it turns out that the dynamical equa- 
tions for the angle field show the proper behavior, that 
is, no long-range interactions between grain boundaries 
arise. Therefore, we believe that this approach can be a 
useful alternative to existing orientation-field models. 

The remainder of the paper is organized as follows. In 
Sec. HU we first outline the model (for isotropic inter- 
faces) and then show how to determine a suitable cou- 
pling function between phase and orientation fields that 
yields localized grain boundaries. Furthermore, we dis- 
cuss in detail the properties of these grain boundary so- 
lutions. In Sec. IIII1 we then extend the model to include 
intcrfacial anisotropy, and present numerical tests both 
in one and two dimensions to show that the model has 
the desired properties. Perspectives for future work are 
discussed in Sec. HVl 



II. MODEL 



General remarks and outline 



We consider a two-dimensional system in which the 
crystal orientation can be described by a single angle field 
#(x, t)\ the local state of the system is described by a 
phase field </)(x, i). We start from a free energy of the 
form 



F 



c/x 



D 



V0) 2 + HV(<}>, T) + Kg(cf>)(Vdf 



(1) 



where D, H et K are constants. Furthermore, the di- 
mensionless function V((j>, T) is a temperature-dependent 
double-well potential with minima at <f> = 1 (solid) and 
4> = (liquid), and g is a function that tends to zero 
for (j> — > 0. For simplicity, we have neglected crystalline 
anisotropy for the moment; its inclusion will be discussed 
below. 

Let us first consider a coupling function g{<f>) that 
takes a finite value g (I) in the solid. The scaling ar- 
gument which proves that there are no localized grain 
boundary solutions in this case is as follows. In the ab- 
sence of anisotropy, the grain boundary solution must 
depend only on the total angle change (misorientation) 
A9 between the two grains. For a solid of size L with 
a homogeneous continuous orientation variation between 
A9/2 and —A9/2, the energy cost (with respect to a 
state of constant 9) is g(l)K L(A9 / L) 2 which decreases 
monotonously with growing L and tends to zero when 
L — > oo. In contrast, any localized grain boundary so- 
lution has an energy that is independent of the system 
size (and includes the energy cost for the variation of 
the phase field). For sufficiently large L, the uniform 
angle variation has therefore always lower energy than a 
localized grain boundary. More precisely, whereas grain 
boundaries can be stable in finite systems depending on 
the values of L and g(l), an isolated grain boundary in 
an infinite system is always unstable. This scaling ar- 
gument also clarifies why a term of the form \V9\ cures 
this problem: with such a term, the energy of a homoge- 
neous variation of the angle is independent of the system 
size and can be made larger than the energy of a local- 
ized solution by a proper choice of coupling functions and 
parameters. 

Here, we propose an alternative approach to achieve 
the same goal. Since the above argument holds only for a 
finite value of <?(!), a possible way out is to use a singular 
coupling function g(<p) that tends to infinity when <p —¥ 1. 
The rationale for this idea is that in a crystalline material 
any continuous variation of the local orientation implies 
an elastic or plastic strain, which has a high energy cost. 
We will show below that we can indeed choose a singular 
coupling function that leads to grain boundaries with the 
desired properties. This function cannot be explicitly 
linked to microscopic physics or an elastic model; our 
approach is therefore essentially phenomenological. 
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For the sake of definiteness, we use a coupling function 
of the form 



g(4>) = 



(i 



(2) 



where a is some positive real exponent and / is a poly- 
nomial in 4> which satisfies /(0) = and /(l) = 1. 

It is useful for the following developments to non- 
dimensionalize the free energy functional and the in- 
volved length scales. The functional has three param- 
eters: D and K have dimensions of energy per unit 
length, and H sets a free energy (density) scale for the 
(dimensionless) double- well potential V(<fi, T). Dividing 
the whole functional by H, we obtain 



F = dx 



y(0,T) + w^ 2 5 (0)(v0) 2 



(3) 

where W$ = \/D/H and We = V 'K/H are the char- 
acteristic length scales associated with solid-liquid in- 
terfaces and grain boundaries, respectively. Scaling all 
lengths with W$ and defining fj, — Wg/W^, we finally 
obtain 



F = I dx 



(V0) 2 + ^,T) +/ x 2 . 9 (0)(V0) 2 



(4) 



In the following, we consider the case \x = 1 that can 
be obtained by a proper rescaling of by /i. It should 
be pointed out that this is also the best case for a use- 
ful model since then the characteristic thickness of solid- 
liquid interfaces and grain boundaries are comparable. 

The equations of motion for the two fields are obtained 
by the standard variational procedure [l3l - [l8j . 



P(cf>, ve)Ttd t <t> 



SF 



6F 



Q(q>,V6)Tgd t e = - — 



(5) 
(6) 



where and rg are relaxation times, and P{4>, V0) and 
Q(<p,V9) are (dimensionless) mobility functions. 

For the solidification of a pure substance, the motion 
of solid-liquid interfaces is linked to heat transport by the 
Stefan condition (the latent heat released or consumed at 
moving interfaces must be transported by the local diffu- 
sion currents). The model therefore has to be completed 
by an equation for the temperature field. We define a 
dimensionless temperature u by 



m(x, t) 



T(x,t) 



L/c 



(7) 



where T m , L, and c are the melting temperature, the 
latent heat of fusion, and the specific heat, respectively. 
The field u obeys the same equation as in the standard 
phase-field model, 



d t u = D th W 2 u + d t (f>, 



(8) 



with D t h the thermal diffusion coefficient. 

It turns out that in order to obtain a model with the 
desired properties, a certain number of conditions on the 
exponent a, the polynomial / and the double-well po- 
tential V have to be satisfied that are linked in a some- 
times subtle way to the internal structure of the grain 
boundary. In order to facilitate the understanding of the 
following developments, it seems useful at this point to 
state a summary of these results, while the details of the 
arguments will be given below. The main points are the 
following. 

1. Localized grain boundary solutions which smoothly 
connect to an infinite bulk solid on both sides can 
exist only for a > 2. 

2. For a fixed total misorientation A9 and increas- 
ing system size L, the energy cost of a uniformly 
strained solid tends to zero for a < 2 and to infinity 
for a > 2. Together with condition 1, this yields 
that the only reasonable value for a is 2. 

3. The requirement that the angle variation is strongly 
concentrated in the center of the grain boundary 
imposes that the first term appearing in the poly- 
nomial / is of order </> 3 . 

4. A uniformly strained solid must be stable against 
the spontaneous formation of grain boundaries. 
This imposes the relation ab — d > 0, where a = 
V"{l)/2, d = -V"'(l)/6, and b = —/'(I), and 
primes denote a derivative with respect to (j>. 

5. For a = 2, the competition between the gradient 
terms of the phase field and the orientation field 
makes it impossible to obtain stable grain bound- 
aries with misorientations smaller than a critical 
misorientation A8 rn . The magnitude of A9 m is con- 
trolled by the same combination of parameters that 
appears in condition 4, ab — d. 

In particular, the last point has implications for the 
choice of the double-well potential V. In the stan- 
dard phase-field model for solidification, this potential 
is V{(j)) = 2 (I - <j)) 2 + uA(IO0 3 - 150 4 + 6<?(> 5 ), where 
A is a dimensionless coupling parameter. For this po- 
tential, the second derivatives with respect to the phase 
field, V"((p,T), are independent of temperature in both 
wells (0 = and <j) — 1); however, the third derivatives 
depend on the temperature. For our model, this would 
yield a temperature-dependent minimal grain boundary 
misorientation. We prefer to use a potential where, in ad- 
dition to the above properties, the third derivative of V 
in the solid well is kept constant. Furthermore, since for 
the calculations of the grain boundary energy the solid 
state is the reference state, we wish to keep the energy of 
the solid potential well independent of temperature and 
equal to zero, whereas the height of the liquid potential 
minimum varies with temperature. 
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A set of model functions which satisfies all the above 
requirements and conditions is 



70 3 - 60 4 



(9) 



-uA(l-2O0 3 + 450 4 -360 5 + lO0 6 ) (10) 



These functions will be used in the numerical examples 
given below. In the following, we will first give the de- 
tails of the arguments that lead to the restrictions on 
these functions, and then discuss their consequences on 
the structure of the grain boundary solutions, in par- 
ticular concerning the dependence of these solutions on 
temperature. 



Construction of the coupling function: the 
details 



Our model should have one-dimensional grain bound- 
aries as equilibrium solutions. In the framework used 
here, a grain boundary is a localized interface (region 
where <fi is smaller than 1) between two solids (semi- 
infinite regions where 0=1) with different crystal ori- 
entations. Therefore, we are looking for appropriate sta- 
tionary solution of Eqs. ([5]) and (J6)). The calculation of 
the functional derivative yields 

= d x (g(^)d x 0), (11) 
= -V'{<f>, T) + d xx <i> - g'{4>){d x 9) 2 . (12) 

From Eq. (jlip. we obtain that at equilibrium 



9 (<t>)d x e = c, 



(13) 



where C is a constant to be determined later. Inserting 
this expression into Eq. (|12p yields 



d xx <j> = V'(<l>,T)+g'(<fi) 



C 2 



c 2 



g(4>) 



(14) 
(15) 



The second equality puts this equation in a form suitable 
for the use of the well-known analogy with a mechanical 
system: replacing x by t and <p by y, one obtains the 
equation of a particle with coordinate y moving in a po- 
tential 



V ef{ (y) = -V(y,T) + 



C 2 

g(y) 



(16) 



(see Fig. [T] (a) ) . The grain boundary solution corresponds 
then to the trajectory of a particle that starts with zero 
velocity at the unstable equilibrium point y = 1, (the 
point A in Fig.QJa), which corresponds to the solid state) 
and reaches the turning point y = yo (the point B in 



Fig. Q] (a)), in an infinite time and comes back. The in- 
finite duration of the trajectory corresponds to the con- 
dition that the limit of 4> in the bulk is <j> = 1 (solid) 
on both sides of the grain boundary. This leads to the 
requirement that </> = 1 must be an equilibrium point of 



off, i-e. 



dV eS 



-V'(l)-C 2 (l 



_ l Ct f + (!-</>) f 
f 2 



where we have supposed a > 1 for V c g to be diffcrcntiablc 
in cj) = 1 . In addition, 0=1 must correspond to an unsta- 
ble equilibrium point of V e s, that is, the second derivative 
of the effective potential must be negative for </> — > 1. A 
simple argument helps to determine the sign of the sec- 
ond derivative. Since the double-well potential V(4>, T) 
has a minimum at (f> = 1, we can expand it around this 
value as V(cj>, T) w V(1,T) + V"(1,T)(1 - </>) 2 /2. Then 
the effective potential becomes 



V* 



-V(1,T)+(1- 



-V"(l,T)/2 + C 



2(1 



\a-2 



(18) 

From this expression, it is clear that the sign of at 
4> = 1 is determined by the sign of the function s{(f) = 
(-V"(l,T)/2+C 2 {l~'(t)) a - 2 /f ■((/>)) at (f> = 1. If at-2 > 0, 
the second term in this function tends to zero for <fi —¥ 
1, whereas — V"(1,T) is finite and negative. Therefore, 
s(l) is negative. In contrast, if a — 2 < 0, the second 
term becomes dominant for <fi — ¥ 1, and s is positive. If 
a = 2, s = (-^"(1, T)/2 + C 2 /f (0)). In this case, since 
f(l) = 1, s changes sign when C 2 crosses —V"(l,T)/2. 

In summary, the fact that 1 must be an unstable equi- 
librium point of V e ff imposes to choose a > 2. We pro- 
ceed by showing that the study of the homogeneous solid 
solution imposes a — 2. To do this, we consider a long 
homogeneous solid of length L where 9 is varying from 
A8/2 to -A9/2 with constant gradient d x 9 = A9/L. 
The scaling argument outlined above shows that a reg- 
ular coupling function yields an energy cost of such a 
homogeneous solution that goes to zero when L — > oo at 
constant A9] we would hence like to construct a model in 
which this is not the case. Qualitatively, it is easy to see 
why this can be achieved with a singular coupling func- 
tion. As long as d x 9 is different from zero, the singularity 
forbids the phase field to take the value <fi — 1 correspond- 
ing to the solid, since this would imply an infinite energy 
cost. However, when d x 9 decreases, the phase field can 
get closer and closer to <fi = 1. Therefore, whereas (d x 9) 2 
decreases, the value of g(4>) increases, which can be used 
to tune the dependence of the total energy on the system 
size. The requirement that the energy must remain finite 
(it should not go to zero, but not tend to infinity either 
when L — > oo) then fixes the exponent a. 

Let us now make this argument quantitative. Since we 
consider a weakly distorted solid, we can assume that cf> 
is close to 1. Let 5 = 1 — 4' where 4> is the (constant) 
value of the phase field. Since (f> is constant, the gradient 
in the phase field does not contribute to the energy of 
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FIG. 1. (a) Schematic drawing of a typical effective potential 
Ves- The starting point of the trajectory corresponding to 
4> = 1 is A, the turning point is B. In both A and B, d x 4> = 
0. (b) Typical profiles of (j> (solid line) and 9 (dashed line) 
through a grain boundary. A corresponds to the limit <j> — > 1 
and x — > ±oo. (c) Grain boundary solution, plotted in the 
phase space spanned by (f> and d x <i>- 



this state. Keeping the dominant term in 5 in the two 
remaining terms of the free energy functional, we obtain 
for the energy E of the state with homogeneous angle 
variation, 



E = L 



(d x ef 



(19) 



Minimization of this expression with respect to 5 yields 

\2\ V(«+2) 



g(d x ey 

V"(l) 

Using this expression in Eq. (|19p , we find 
E L Q (-) 



(20) 



(21) 



where Q = [V"(l)/2][a/V"(l)} 2 ^ a+2 '> + [V"(l)/a} a ^ a + 2 '> 
is positive. The only possible value for a that yields a 
finite energy in the limit of large L is a = 2. With this 
choice, the energy of a large system with constant angle 
variation is proportional to AO. 

The requirement that this homogeneous solution must 
be stable yields another condition. To understand its 
origin, consider the curve E(L). If this curve is concave, 
there is an instability that is analogous to spinodal de- 
composition: at fixed total AO, it is more advantageous 
for the system to concentrate the distortion in one or sev- 
eral regions because this lowers the total energy; in other 
words, there would be spontaneous formation of grain 
boundaries. To prevent this, the curve E(L) should be 
convex at least for low values of the distortion (weakly 
strained solids), which implies that it has to tend to the 
limit value from above. To find the corresponding math- 
ematical condition, the expansion in 5 has to be pushed 
to higher order, where the energy writes (for a = 2) 



E = L 



a5 2 + dS 3 + — 



AO 
L 



1 + bS + cS 2 



(22) 



where a = V"(l)/2, d = -V"'(l)/6, b = 
c = /"(l)/2. Using the fact that {A6/L) 2 



/'(I), and 
- S a+2 ac- 
cording to Eq. pn)l . it is straightforward to show that 
the requirement of E(L) to be a decreasing function of L 
for L — > oo implies 



ab + d > 0. 



(23) 



To finish with the study of the homogeneous solution, 
we find it useful to present some numerical results using 
the complete model. The computation of the energy as 
a function of L for a fixed value of AO gives the expected 
result that the energy density depends only on the value 
of A0/L : E/L = £(A0/L). The function £(A0/L) is 
plotted in figure [2] (a). It is interesting to note that in 
this case, the value of dE/dL is a function of A9/L only. 
These results also show that, up to a rescaling of the 
energy E and the length L, all the curves of the energy 
as a function of L collapse onto a single master curve 
which is plotted in figure [2] (b). This curve presents a 
maximum, the position of which is given by the solution 
of the equation 

£(A0/L) = £'{A0/L). (24) 



C. Grain boundary solutions 

We now turn to a more complete description of the 
grain boundary equilibrium solutions for a > 2. They 
can be easily computed using the mechanical analog out- 
lined above. For a given value of C, one first computes 
the value of 4>o corresponding to the turning point (B in 
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FIG. 2. (a) solid line: E/L = £(A8) as a function of 
V# = A8/L, dashed line: (j> as a function of A9/L. For 
values of |V6>| > 0.119, the solid state with constant orienta- 
tion gradient ceases to exist, and the only possible solution 
is the liquid state. The undercooling is taken equal to 0.1 
(it = -0.1). (b): curve E = f(L) for A8 = .1 for the same 
parameters as in (a). Inset: blowup of the vicinity of the 
maximum. 



fig. [T]) by solving V c g(l) = V c s{4>o)- Then, the equilib- 
rium solution (going from the middle of the grain bound- 
ary to the solid state on one side of the grain boundary) 
corresponds to the trajectory of a moving particle in V c g 
that starts at the point B with zero velocity. Once 4>{x) 
for the equilibrium solution has been computed, 6{x) is 
given by Eq. (jlip . and the total angle variation through 
the grain boundary can be obtained by simple integra- 
tion. In order to compute the grain boundary solution for 
a given misorientation, the proper value of C has to be 
found, for instance by using a Newton-Raphson method. 

In figure [3Ja) we present a grain boundary solution for 
the same total angle variation A9 — 0.6 for a = 2 and 
a = 3 (in the calculations for a = 3, we use the cou- 
pling function g of Eq. @, changing just the exponent). 
As predicted, the angle variation is concentrated in the 
vicinity of the minimum of <j) (the turning point B), and is 
very small outside the region where <\> differs from 1 . This 
indicates that this is indeed a localized grain boundary 
solution. 

Using the computed solutions, the grain boundary en- 
ergy is given by 



E, 



GB 



d.r 



(V</>) 2 + V(^T)+g(cp){V0f 



(25) 
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FIG. 3. (a): typical grain boundary profiles for the same mis- 
orientation (A8 = 0.6) and two values of a: a = 2 (full lines) 
and a = 3 (dashed lines). The dimensionless undercooling is 
0.1. (b): Grain boundary energy as a function of the misori- 
entation A6 for different values of a (a = 2, 3 and 4). The 
inset displays a zoom around the origin and shows the sin- 
gular behavior for a = 3 and 4 and the fact that for a = 2 
there is a small band of misorientations for which no grain 
boundary solution exists. (Here fi 2 = 10) 



out any angle variation as the zero of the free energy 
density). We display this energy as a function of the to- 
tal angle variation AO in Fig. [3] (b). In the cases a = 
3 and 4, for small values of A6, the energy has a power- 
law behavior (E GB - (A6») 7 , with 7 =0.8 and 0.66, re- 
spectively). When AO increases toward infinity (this is 
not absurd since we have previously rescaled AO by /i), 
the grain boundary energy converges toward a constant 
value. 

In the case a = 2, the grain boundary energy also con- 
verges toward a constant value when AO becomes large. 
However, the behavior at small misorientations is very 
different: for small values of AO, the grain boundary en- 
ergy decreases linearly with the misorientation, but there 
exists a minimal value A0 m for which a grain boundary 
solution can exist (see inset). In fact, in the following we 
show that there is no grain boundary solution with angle 
variation smaller than this value in the case a — 2. 

As mentioned earlier, for a = 2, <fi = 1 is not a max- 
imum of the effective potential V e g for all values of C. 
More precisely, V e g can be formally written as 



(recall that we have chosen the homogeneous solid with- 



KffW>) = (1 - 0) 2 [-p(^, T) + C 2 //(0)], (26) 
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where p(T,4>) is a function such that p(<j>, T)(l — </>) 2 = 
V((j),T). Using the notations of Eq. ([22]) . V c s writes in 
the vicinity of <fi = 1: 

^)=^ 2 (-a-^ +IT ^) (27) 

with (5=1 — 0. It can be seen that the curvature of 
V e f{ changes sign in <j> = 1 when C becomes bigger than 
a. This means that <fi = 1 is now a local minimum of 
V e fi- This change can happen in two ways: either the 
effective potential becomes greater that zero for all <fi < 1 
(the turning point B in Fig. [TJa) moves to the right of 
the point A), or a new maximum develops between the 
points A and B (this happens in our case, as illustrated 
in Fig. @|, which means that a zero of V e s must cross the 
point A from above. It is straightforward to show that 
the latter always happens if ab + d > 0. This explains the 
threshold value for A8: grain boundary solutions with 
finite A.8 exist until the new maximum appears. Beyond 
this point, no grain boundary solution connecting to <f> = 
1 is possible any more. 

We have also verified that all the grain boundary so- 
lutions found up to now can be reached by simple inte- 
gration of the equations of motion, Eqs. ([5]) and ([6]), at 
constant temperature (the field u is held constant), when 
the system is started from a homogeneous solid (with a 
value of the phase field slightly lower than unity) and a 
step function in the angle field. The resulting station- 
ary solutions as well as their energies are independent of 
the system size (for sufficiently large systems) and of the 
discretization used (for grid spacings Ax smaller than 
about 0.5 W). This shows that the model formulation 
is robust, and that the basin of attraction of the grain 
boundary solutions is large enough to make the model 
suitable for simulations of polycrystals. 

D. Temperature dependence of the grain boundary 
solutions 

In order to better characterize the model, we find it 
useful to explore the temperature dependence of the grain 
boundary solutions. We start out with the regime that 
should be of most practical interest, namely temperatures 
below the melting point (T < T m , u < 0). 

From the discussion of the preceding section, it is clear 
that the minimal grain boundary angle strongly depends 
on the behavior of V(<j>) close to the solid state. Since we 
have chosen V(<j>) such that the solid well is independent 
of temperature, the small-angle solutions should not be 
much affected by temperature changes. The same argu- 
ment indicates that changes in undercooling (increasing 
undercooling) should only appreciably modify the high- 
angle grain boundaries and their energies. Simulation 
results show that this is indeed the case (see Fig. [5]). 
As expected, for the double-well function of Eq. (fTTTj) 
that was designed to maintain the combination ab — d 
constant (see appendix A), the minimal grain boundary 



0.94 0.95 




FIG. 4. Effective potential V e a for a = 2 and two values of C: 
C = 1 ± 5 x 10~ 4 (dashed and solid lines, respectively). The 
thin lines correspond to the interval 0.94 < (j> < 0.95 (upper 
axis) and display the turning point of the effective potential. 
The thick lines correspond to the interval 0.999 < <j> < 1.001 
(lower axis) and display the vicinity of the solid state. The 
arrows indicate the direction of change when C is lowered 
It can clearly be seen that the turning point undergoes no 
qualitative change, whereas a new maximum appears close to 
4> — 1 due to the fact that a zero of V e s has crossed the point 

4> = i. 



angle A9 m is found to be independent of undercooling. 
The grain boundary energy for small angles is also al- 
most independent of temperature (see Fig. [5] (b)). For 
high-angle grain boundaries, the grain boundary energy 
depends roughly linearly on the undercooling, which is 
expected since for such boundaries the center of the grain 
boundary is almost in the liquid state (</> « 0), and the 
free energy difference between solid and liquid is propor- 
tional to T-T m . 

We now turn to the case of negative undercooling 
(overheated solids). Whereas it is difficult if not impos- 
sible to overheat solids in experiments, due to the fun- 
damental assymetry between solidification and melting 
(for a review, see 24]), the existence of grain bound- 
ary solutions above the melting point has recently been 
actively discussed in the context of grain-boundary pre- 
melting, due to new results coming from molecular dy- 
namics [25|, [26| and phase-field crystal [13, HI] simula- 
tions. It has also been examined in orientation-field 
[i~5l [l|| and multi-phase- field models [29M3l| . It seems 
therefore useful to analyze our model in the light of this 
discussion. 

A grain boundary can formally be described as two 
solid-liquid interfaces separated by a thin disordered (al- 
most liquid) layer. The behaviour of a grain boundary 
when the melting point is approached then depends on 
whether the grain boundary energy is larger or smaller 
than the excess free energy of two solid-liquid interfaces. 
In the former case, the grain boundary is repulsive, since 
at the melting point, where the free energy density of 
both phases is equal, it is advantageous to separate the 
two solids by a macroscopic liquid layer. In the latter 
case, the grain boundary is attractive, since it is favor- 
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able for the two solids to stick together. More quantita- 
tively, the total free energy excess of the grain boundary 
can be described as twice the solid-liquid surface free en- 
ergy plus a so-called disjoining potential, which describes 
the interaction between the two interfaces as a function 
of their distance. In the classic theory of grain bound- 
ary premelting, this disjoining potential is assumed to 
be a simple exponential (see for example |29j), but it 
was shown that this simple form is not sufficient to de- 
scribe low- ang le grain boundaries in the phase-field crys- 
tal model [28|]. Indeed, for these boundaries the interac- 
tion is attractive for large separations, but repulsive for 
short ones. As a consequence, equilibrium grain bound- 
aries of finite width can exist at and above the melting 
point. It was later shown that such behavior (long-range 
attraction and short-range repulsion) is als o g eneric for 
the standard multi-phase-field models (see [31| for a de- 
tailed discussion). 

In our model, the structure of the solution space can be 
directly and simply deduced from the effective potential 
and the relation between the constant C and the total 
angle variation A9. Let us start out by a brief discus- 
sion of the solutions that exist when the angle variation 
is zero. In this situation, d x 9 is zero everywhere, which 
implies C = 0, and the effective potential is just the 
negative of the double- well potential. For T < T m , the 
solid has a lower free energy than the liquid, which im- 
plies that the solid maximum in the effective potential is 
higher than the liquid one. Therefore, no turning point 
exists, and no grain boundary solution of the type out- 
lined above is possible. For T > T m , the liquid maximum 
in V e s is higher than the solid one, and therefore there 
always exists a turning point. The corresponding grain 
boundary solution is unstable. Indeed, it can be seen as 
a thin layer of liquid sandwiched between two solids, in 
a situation where the attraction between the two solid- 
liquid interfaces is exactly balanced by the free energy 
difference between solid and liquid. Hence, this solution 
corresponds to a saddle point of the free energy, which 
can be seen as a one-dimensional liquid nucleus. When 
the melting point is approached from above, the thick- 
ness of the liquid layer diverges, and the grain boundary 
energy approaches twice the solid-liquid surface tension. 

Next, consider the situation when T < T m . If the an- 
gle variation is finite, the constant C differs from zero, 
and the effective potential has two contributions. Since 
.9^0 for <j) — > 0, the second term in Eq. (Tl6|) di- 
verges in this limit. In addition, since g diverges for 
4> — >■ I, l/g vanishes for <j> — »• I. Since C 2 /g is pos- 
itive, there is necessarily at least one solution for the 
equation V c s(<f>o) = Vtff (1), which implies that there is a 
turning point for any value of C . It remains to link the 
values of C to the total angle variation. Whereas we did 
not find an analytical expression for this relationship, the 
main trends can be extracted from an analysis of the tra- 
jectory close to the turning point. When T tends to the 
melting temperature, the constant C necessarily becomes 
small since the liquid maximum in the returned double- 



well potential lies closely below zero. This implies that 
the turning point is very close to the liquid state. As a 
consequence, g is also small, and the orientation gradient 
becomes large. Since g behaves like a power law with 
exponent larger that one close to the liquid state (see 
appendix A), the gradient at the center of the boundary 
becomes larger when C becomes smaller. Therefore, ar- 
bitrarily large total angle variations can be achieved for 
any temperature below T m . 

This situation changes when T > T m : now, the liq- 
uid maximum in the returned double-well potential is 
higher than zero, which means that a turning point ex- 
ists even for C = 0. When a non-zero value is chosen 
for C, this turning point moves further away from the 
liquid state. Therefore, the phase field remains finite in 
the center of the grain boundary, and the orientation gra- 
dient in the center of the grain boundary tends to zero 
when C — > 0. In this case, the curve of A9 versus C is 
non-monotonous and presents a maximum. This implies 
that there is a maximum misorientation beyond which no 
grain boundary solutions exist. This maximum misorien- 
tation depends on the temperature and tends to infinity 
when T — > T m . For misorientations smaller than this 
critical value, there are two distinct solutions. The so- 
lutions for C smaller than the maximum correspond to 
unstable liquid nuclei, whereas the solution branch for C 
larger than the maximum connects to the stable grain 
boundary solutions discussed previously. Solutions cor- 
responding to a grain boundary and an unstable liquid 
nucleus are shown on Fig. [5] (d) . While the angle profile 
is extremely similar in both cases, the liquid nucleus so- 
lution has a much smaller minimal <f>. When the energies 
are plotted versus misorientation, the two branches meet 
at a cusp, as shown in Fig. [5] (b). At a fixed misorienta- 
tion above the minimal value discussed previously A9 m , 
there are thus stable grain boundary solutions for any 
T < T m , and for a range of temperatures T > T m that 
depends on the misorientation. This overheated range 
tends to zero when the misorientation becomes large. 

The stable grain boundary solutions all exhibit a devi- 
ation of the phase field from the solid state (<j) = I), which 
can be interpreted as a thin layer of liquid by performing 
a Gibbs construction as in Ref. [28| . The existence of such 
a layer of finite thickness results from the interplay be- 
tween two antagonistic effects, which can be understood 
through the interpretation of the grain boundary as two 
solid-liquid interfaces that come close together. On the 
one hand, the overlap of the two phase-field "tails" of 
the two interfaces creates an attractive force that would 
eliminate the liquid layer (as is happens in the case of 
zero misorientation). On the other hand, the finite total 
misorientation generates a restoring force that opposes 
this compression of the liquid layer: if the two interfaces 
get closer together, the orientation gradients have to in- 
crease. On the branch on solutions that is stable, this 
increase of the gradient leads to an increase of the free 
energy of the grain boundary solution, whereas the op- 
posite is true on the unstable branch. 
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FIG. 5. (a) : Grain boundary energy as a function of A8 for 
different values of the undercooling. The dash-dotted curve 
corresponds to u — 0. (b) : Grain boundary energy as a func- 
tion of misorientation for different values of the undercooling 
for overheated grain boundaries (u > 0, negative undercool- 
ing) . Grain boundary solutions exist only for a finite range of 
misorientations. (c) : minimal value of <j> in the grain bound- 
ary in the case of overheated grain boundaries, as a function 
of AO, inset: value of the maximal angle variation through a 
stable grain boundary as a function of tt -1//4 . (d) : (p (solid) 
and 8 (dashed) as a function of x along both a grain bound- 
ary and a one-dimensional liquid nucleus for the same angle 
variation A9 = 0.4 and u = 0.005. 
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FIG. 6. Grey scale plot of <j> for an equilibrium trijunction. 
The three grains with respective grain orientation #0, —80 and 
correspond to the white regions, and the grain boundaries 
correspond to the dark region. Here, 80 = 10° , the undercool- 
ing was 0.085, and fi 2 = 11.2. 



with an angle variation bigger than a threshold value 
(that can be made as small as desired by a proper choice 
of the quantity ab + d). It also has solutions that cor- 
respond to a homogeneous solid with constant gradient 
in crystal orientation as long as this gradient is smaller 
than a threshold value. For a given total angle variation, 
the free energy of the grain boundary is always smaller 
than the energy of the homogeneous solution (which may 
be metastable). 



As a result, all the grain boundaries in our model 
are attractive. This finding coincides with the classic 
criterion. Indeed, the solid-liquid surface tension can 
be calculated analytically in our model and is equal to 
7sl = a/2/6. The curve of Eqb versus AO at u = 
(T = T m ), shown in Fig. El has an asymptote that tends 
exactly to 27 s i = v2/3 ~ 0.4714 from below, which cor- 
responds to two solid-liquid interfaces with a weak at- 
tractive interaction that is mediated by the exponential 
tails of the interfaces (see below). This behavior cannot 
easily be changed in our model. For example, the method 
used in Ref. [3l| to make the interfaces repulsive is not 
applicable here because it requires a free energy density 
with minima in both field variables, which is not allowed 
if 6 is to be interpreted as an orientation and rotational 
invariance is desired. 

More complex behavior has been found in other 
orientation-field models [HI, [l8| by a suitable tuning of 
the various coupling functions. In particular, for certain 
temperature ranges the possibility of coexistence between 
two different grain boundary states was found. Similar 
behavior could probably be generated in our model if the 
double-well potential and the coupling function g(cj>) are 
modified, but this issue is not pursued any further here. 

To summarize, the model presented here with an ap- 
propriate choice of the coupling function and double- well 
potential has stable localized grain-boundary solutions 



E. Trijunction points and Young's law 

We have also checked explicitly that Young's law is 
satisfied at trijunction points. In order to obtain a static 
trijunction and to avoid curvature effects induced by sim- 
ple Neumann boundary conditions, the simulation was 
performed as follows: first, we let three initially circular 
grains grow in a simulation box with Neumann bound- 
ary conditions. When the grains have impinged on the 
system boundaries, we use the value of the fields at the 
boundary at a certain time as Dirichlet boundary condi- 
tion for the further evolution, which imposes that neither 
the phase field nor 6 evolves any further at the bound- 
ary. This amounts to pinning the grain boundaries at 
fixed points of the system boundary. For initial orien- 
tations of the three grains given by 0q, 0, and — 6>o as 
displayed in Fig. O Young's law predicts 

{ E GB (29 ) \ 

Vj = acos (^m) ' ^ 

where ip is the half of the dihedral angle, as sketched in 
Fig. |H1 Egb(&o) and Egb{2&o) are the energies of grain 
boundaries with misorientations 9q and 28q that have 
been calculated in the previous subsection. For #0 = 
4.2° (resp. 10°), the angle ip was 46° (resp. 55°), to be 
compared to 46° (resp. 52.4°) predicted from Young's 
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FIG. 7. Top: Sketch of the tricristal configuration, which 
can evolve towards a homogeneous system by the rotation of 
the central slab. Bottom: Rotation rate of the central slab 
as a function of distance between the boundaries for two dif- 
ferent values of the angle difference through the grain bound- 
ary (AO) (left) and as a function of misorientation (AO) for 
a given value of the distance between the two grain bound- 
aries (right). In these simulations, the orientation of the outer 
grains is kept constant. 



law. The agreement is excellent, which is not surprising 
because Young's law should generally be valid in phase- 
field models that derive from an energy minimization [32j . 



F. Grain boundary tails and the tricristal 

Next, we briefly examine the interaction of two grain 
boundaries in the tricristal configuration sketched in 
Fig. [7] two grains of the same orientation sandwich a 
slab of material with a different orientation. The two 
grain boundaries obviously contribute a finite free energy. 
Since the evolution of the system is variational according 
to Eqs. ([3]) and ^ and nothing in the equations for- 
bids the central crystal slab to rotate, the system tends 
to eliminate these defects: the orientation of the central 
slab slowly changes with time to approach the one of the 
outer grains, until the angle field finally becomes homoge- 
neous. Obviously, for the sketched geometry this change 
cannot correspond to rigid body rotation. The evolution 
equation for the orientation actually describes the local 
rotation of disconnected objects, such as the molecules of 
a liquid crystal. For a continuously connected rigid solid, 
such a rotation could only be accomplished by a change in 
local connectivity, obtained for example by the migration 
of dislocations from one grain boundary to the other (see 
[33j for a more detailed discussion). However, since this 
process (if it takes place at all) would be extremely slow, 
the rotation of the central slab should be suppressed. In 
the previous orientation-field models that contain a mod- 



ulus of the gradient in the free energy functional, this is 
problematic since the functional derivative has a singu- 
larity at the points with vanishing orientation gradients 
- that is, inside the grains. The singular diffusion equa- 
tion that results from the functional derivative leads to 
a non-local coupling of the grain boundaries 22]. As a 
consequence, the rate of rotation of the central crystal 
slab does not depend on the distance between the grain 
boundaries. Therefore, the only solution to suppress the 
unwanted rotation is to choose a mobility function which 
vanishes inside the crystal grains (lfij . 

In our model, the interaction between the grain bound- 
aries is local. This can be easily demonstrated by an 
analysis of the "tails" of the boundary, that is, the ap- 
proach of the phase and orientation fields to their bulk 
values. This is again straightforward by using the me- 
chanical analog. Indeed, a linearization of Eq. (|15p in 
the vicinity of the solid state yields 



d 2 V cS 



(0-1), 



the solution of which is 



1 — Aexp[±^(x — xq)] 



(29) 



(30) 



and the values of the in- 



with£ = ^(d 2 K ff )/(d^)| 0=1 , 

tegration constants A and xq as well as the sign inside the 
exponential are fixed by the boundary conditions (which 
depend on the position of the grain boundary) . The gra- 
dient in orientation is then obtained from Eq. (1131) and 
tends to zero in the grains since g(4>) diverges as cf) — >■ 1. 

Several conclusions can be drawn from this result. 
First, the grain boundaries in our model are truly lo- 
calized, in the sense that in a heterogeneous system only 
an exponentially small part of the orientation variation 
takes place outside the grain boundaries. Second, the in- 
teraction between two grain boundaries decays exponen- 
tially with their distance for well-separated grain bound- 
aries, since for large separation the interaction is medi- 
ated by the overlap of the exponential tails (see 3l| for 
a more detailed discussion). As a consequence, even for 
a constant mobility in the orientation evolution equation 
(Q(<fi,\79) — 1), the rotation rate of the central slab in 
the tricristal configuration should rapidly drop with the 
size of the slab. Finally, the characteristic length of decay 
for the exponential tails l/£ is the inverse of the second 
derivative of the effective potential, taken in the solid 
state. As discussed above, the second derivative depends 
of the value of the constant C and vanishes at a critical 
value of C that corresponds to the minimal misorien- 
tation A9 m . Therefore, grain boundaries with smaller 
misorientations interact more strongly. 

These predictions are borne out by simulations on the 
tricristal configuration, as shown in Fig. [7] For grain 
boundaries of a fixed misorientation, the rotation rate of 
the central slab decays exponentially with systems size 
(and thus, separation of the grain boundaries), and the 
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rotation rate is larger for smaller misorientations at a 
fixed separation. Therefore, while our model does per- 
mit grain rotation, qualitatively by the same mechanism 
than in previous orientation-field models, the artefacts 
introduced by this design are small even for a constant 
mobility function, and it can be hoped that they will not 
drastically affect the results of large-scale simulations. 



III. ANISOTROPIC INTERFACES 

Our model can describe two different types of bound- 
aries: solid-liquid interfaces and grain boundaries, both 
of which are in general anisotropic. For solid-liquid inter- 
faces, the surface tension and interface mobility depend 
on the orientation of the interface with respect to the 
crystal axes of the solid. For grain boundaries, the energy 
and mobility depend on the orientation of both crystals 
that meet at the boundary. This yields two independent 
parameters in two dimensions, and five independent pa- 
rameters in three dimensions. Grain boundary energies 
and mobilities are conventionally given as functions of 
the mis orientation between the two crystals and the in- 
clination of the grain boundaries. 

The standard way of including the anisotropy of the 
solid-liquid interfaces in phase-field models is to make 
the coefficient of the phase-field gradient term dependent 
on the interface orientation (see for example [3J|)- We 
will now explore the effects of introducing this type of 
anisotropy in our model. As we will see below, the modi- 
fication of the phase-field gradient term also makes grain 
boundaries anisotropic. Since the anisotropies in both 
types of boundaries are created by one single term in the 
free energy functional, they are obviously not indepen- 
dent. A possible way to choose both anisotropies inde- 
pendently would be to introduce a suitable dependence 
on orientation in the coefficient of the orientation-field 
gradient term. However, this possibility will not be ex- 
plored further here. 

For a crystal of a cubic material in two dimensions with 
its crystallographic axes aligned with the coordinate sys- 
tem, the orientation-dependent surface tension is usually 
written as 

j{<p) =7d [1 + 64(508(4^)], (31) 

where (f is the angle between the interface normal (point- 
ing into the liquid) and the x axis. If the crystal is rotated 
by an angle this expression becomes 

7(^0)= 7 o[l + e 4 cos[4(^-0)]]. (32) 

To implement this anisotropy in the phase-field model, 
we define the unit normal vector by n = — V0/|V(/>| and 
choose the gradient coefficient to be 

W4n,6) = W a(n,6) (33) 

with the anisotropy function 

a(n, 6) = 1 + e 4 [cos(40) (4(n* + n\) - 3) 

+ sin(4#) (An x nl - An y n 3 x )] , (34) 



where we have repeatedly used the addition theorems for 
the trigonometric functions. Lengths are now scaled by 
the average interface thickness Wo ■ 

A few consequences of the fourfold symmetry should 
be mentioned here. The fundamental issue is that for 
a cubic material all orientations that differ by integer 
multiples of tt/2 are equivalent. Indeed, a crystal unit 
cell can be locally rotated by a multiple of tt/2 without 
changing the state of the crystal. Therefore, 9 is actu- 
ally defined modulo tt/2, which means that the values of 
6 can be restricted to the interval [0,7r/2[. This means 
that, contrary to what we did in the previous section, we 
cannot rescale 9 to make fi = 1. Nevertheless, we will 
consider in the following the case \x = 1 for simplicity. 
Furthermore, care has to be taken in the calculation of 
differential operators. For instance, in the evaluation of 
the gradient energy and in the Laplace operator, differ- 
ences between the values of 9 at neighboring grid points 
need to be computed. We follow the method proposed in 
Ref. [2(j: in the evaluation of such operators for a given 
grid point i, we check for each neighboring point which of 
all the equivalent values of 9 gives the lowest value of the 
gradient energy, and calculate the evolution of point i us- 
ing this value. This procedure ensures that local "jumps" 
of the angle field with amplitude tt/2 (or multiples of it) 
do not introduce any energy cost. Thus, all the crys- 
tal symmetries are correctly implemented. Nevertheless, 
for simplicity of exposition, we will use the standard no- 
tation for differential operators in the remainder of the 
paper. 

The evaluation of the variational derivatives in Eqs. 
© and ([6]) yields 

P(^V9)d t <P = V- [a(n,9) 2 V<f\ 

da(n,9)' 



+ d x |V0| 2 a(n,/ 







2 a(n, 



d{d x cj>) 
da(n,9) 



0(M 

V'^,T)-^{Vefg'^) (35) 



and 



Q{<t>,ve)d t 



^V.( g ^)V9)~a(n,9)^^l(V^. 

(36) 

For a uniform (constant) orientation field, the first of 
these equations is identical to the standard anisotropic 
phase- field equation of Ref. [34[ , with the crystalline axes 
rotated by the angle 9 away from the x axis. The equa- 
tion for the orientation field contains the derivative with 
respect to the angle of the phase- field gradient coefficient. 
Physically, this term represents a "torque" exerted on the 
orientation field at the surface. Indeed, if the solid-liquid 
interface is anisotropic, this term drives the angle field 
towards an orientation which corresponds to a minimum 
of the surface free energy. 

In the following, we will examine in more detail the 
influence of the new terms on grain boundaries, solid- 
liquid interfaces, isolated monocrystals and polycrystals. 
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A. Grain boundaries 

We have systematically calculated stationary one- 
dimensional grain boundary solutions of the anisotropic 
equations (|35|) and (|36|) . and determined the grain bound- 
ary energy. The results are presented in Fig. |SJ As 
stated above, in two dimensions, there are two indepen- 
dent parameters. For two grains of respective orienta- 
tions 9 1 and 6* 2 , we define the misorientation, as before, as 
A9 = 92~0i, and the median angle as 9 mi a = (9 1 + 9 2 )/2. 

The curve of the grain boundary energy versus the me- 
dian orientation for a fixed misorientation exhibits the 
expected behaviour, that is, a 4(9 dependence shown by 
the 7r/2 periodicity of the grain boundary energy as a 
function of 6* mi( i. We limit ourselves to inclinations that 
are smaller than 7r/2 since other values can be reached by 
symmetry. It is also interesting to note that higher mis- 
orientations can have a lower grain boundary energy for 
certain inclinations. This is due to the fact that a grain 
boundary consists, roughly speaking, of two solid/liquid 
interfaces and that changing the inclination rotates both 
of these interfaces with respect to the orientation of the 
crystalline axes. This can lead to decrease of the ener- 
getic cost that can be higher than the increase due to an 
increase of the misorientation. 

It should be mentioned that in our first implementa- 
tion of the model, we noticed a very slow motion of one- 
dimensional anisotropic grain boundaries (35j | . This is 
clearly unphysical since the free energy density in the 
grains on both sides is identical, which precludes a per- 
sistent motion to one side. We found that this effect was 
due to the discretization. Indeed, the anisotropic terms 
lead to assymetric <f> profiles at the grain boundary since 
the grain orientation relative to the grain boundary nor- 
mal is, in general, not the same on both sides of the 
interface. Therefore, the unavoidable discretization er- 
rors present on both sides of the boundary do in general 
not exactly compensate each other, which then leads to 
an unwanted movement of the interface. We eliminated 
this problem by changing the discretization scheme: in- 
stead of discretizing the equations of motion (1351) and 
(f36|) . we discretized the free energy functional and re- 
placed the functional derivatives by ordinary derivatives. 
As a result, we obtain discrete equations that derive from 
a discrete energy function. It can easily be shown (see 
appendix B) that with this scheme the energy function 
is strictly decreasing with time; indeed, no spurious mo- 
tion of the grain boundaries took place with this modified 
scheme. 



B. Solid-liquid interfaces and equilibrium crystal 
shapes 

Next, we turn to the study of solid-liquid interfaces. 
Bulk liquid and solid coexist at the melting point (T = 
T m , u = 0), which implies that there must exist an equi- 
librium solution for a planar interface. For a constant 




J mid 



FIG. 8. Top: Grain boundary energy as a function of the mis- 
orientation for different values of the median angle (e = 0.05, 
u = —0.1. Bottom: grain boundary energy as a function of 
the median angle for two different values of the misorienta- 
tion. Inset: same curves, normalized by the minimal energy 
of such a GB (« 0.5 for the case where the misorientation is 
1.5° and ~ 3.2 for a 45° misorientation). Clearly, for small 
values of the misorientation, the grain boundary energy is less 
sensitive to anisotropy than for high values. 



orientation field (0(x) = 9o), this solution (for an inter- 
face normal to the x direction and centered at position 
xq, with the solid located in the domain x < xq) can be 
obtained analytically and reads 



1 — tanh 



(V2W*(06), 



X — Xq 



(37) 



This interface has a surface free energy of 7 — 
'YoW$(6o)/Wa, which has already been used to choose 
Wcf, in Eq. ©. 

However, a constant orientation field is a stationary 
solution of the anisotropic orientation equation (1361) only 
if the angle 9q corresponds to symmetry direction of the 
crystal (a minimum or maximum of the anisotropy func- 
tion). In all other cases, the "torque" term is non-zero 
inside the interface, which generates an evolution of the 
orientation if the system is started from a constant ori- 
entation field. 

This effect is actually not physical for a crystalline ma- 
terial. It is completely analogous to the anchoring ef- 
fect known for liquid crystals. Indeed, in such materials, 
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the constituent molecules are anisotropic and exhibit an 
orientational order, but are free to rotate individually. 
Therefore, if the surface energy depends on the orienta- 
tion of the molecules with respect to the surface normal, a 
competition between bulk and surface effect takes place. 
Any variation of the orientation in the bulk is energet- 
ically penalized, but the energy gain due to the surface 
terms can offset this penalty for molecules close to the 
surface, such that the orientation tends to an energeti- 
cally favorable direction close to the surface. In contrast, 
in a crystalline material, the structure cannot change by 
rotation of the individual crystalline unit cells. The in- 
terface orientation can change only by attachment and 
detachment of atoms. 

To eliminate this effect, existing orientation-field mod- 
els adopt a radical strategy [17|: they just set the 
"torque" term to zero and use Eq. (|36|) with only its first 
term on the right hand side. Whereas this procedure can 
be justified by the arguments given above, it breaks the 
variational structure of the model, which means that sev- 
eral desirable mathematical properties are lost (in partic- 
ular, the free energy is not guaranteed to decrease with 
time any more). Instead of adopting the same strategy 
in our model, we will proceed by showing that the effects 
induced by this surface term are actually small. Once 
again, the facts that the diffusion part of the orientation 
equation is regular and that the orientation variation is 
localized in the vicinity of the interfaces are crucial for 
this property of the model. 

When we start the two coupled equations (l3"5j) and (j3f))) 
from the solution given by Eq. (|37[) at a fixed temperature 
u = 0, we observe that the solid grows very slowly at the 
expense of the liquid. This can be qualitatively under- 
stood by the following reasoning. As mentioned above, 
the torque term induces an evolution of the orientation 
field towards a more favorable direction. If the solid could 
"turn" by changing its orientation, the system would be 
able to reach its global energy minimum (the orientation 
field aligned with a direction of minimal surface tension) . 
However, since the gradient of orientation tends exponen- 
tially to zero away from the interfaces, there is actually 
no term in the orientation equation which could induce 
such a rotation. However, the system can lower its energy 
by moving the interface and simultaneously changing the 
orientation field inside the interface. In this way, a ho- 
mogeneous liquid is successively replaced by a solid with 
a very small "frozen in" orientation gradient that results 
from the time evolution of the surface orientation. This 
"strained" solid has a higher free energy density than the 
liquid, but this can be offset by the lowering of the sur- 
face free energy, such that the total energy of the system 
still decreases. 

In summary, the anisotropy of the solid-liquid interface 
induces an unphysical growth of the solid at the melting 
temperature. In other words, the presence of the ad- 
ditional degree of freedom (change in orientation) leads 
to a shift of the melting point for anisotropic interfaces. 
To quantify the importance of this effect, we have deter- 
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FIG. 9. (Color online) The equilibrium shape of an 
anisotropic crystal (green circles) agrees well with the Wulff 
shape (inner envelope of the red lines). The anisotropy is 
e 4 = 0.05. 



mined the modified melting temperature by performing 
simulations at constant volume of solid, as described in 
appendix C. The shift is actually very small, of the or- 
der of u ~ 10~ 4 . Moreover, since the profile of phase 
field and orientation across the interface are not exactly 
given by Eq. {37}, the actual surface tension slightly dif- 
fers from the intended value. To test the magnitude of 
this effect, we have determined the equilibrium shape of 
a two-dimensional anisotropic crystal, again by perform- 
ing simulations at constant volume of solid. In Fig. El we 
show the resulting equilibrium shape and the compari- 
son to the Wulff construction, which is very satisfactory. 
Indeed, for typical undercoolings that can be reached in 
simulations (u ~ 0.1), the shift of the equilibrium point 
is negligibly small, and therefore it is not surprising that 
the model is capable of producing an excellent approxi- 
mation to the equilibrium crystal shape. 

We have not explicitly checked Young's law at trijunc- 
tion points for the anisotropic version of our model. How- 
ever, it is shown in Ref. [32| that Young's law is generally 
valid for phase-field models that are variational, since 
equilibrium states result from an energy minimization. 
Since this is the case for our model, and the equilibrium 
crystal shape is well reproduced, we are fairly confident 
that our model also correctly implements the balance of 
Herring torque terms at trij unction points. 
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FIG. 10. Snapshots of the numerical simulation of the 
isothermal growth of multiple grains. The grey scale indicates 
the orientation of grains. The undercooling is kept homoge- 
neous and constant and the initial grains (top left) grow until 
they begin to interact with each other (top right). Once all 
the liquid is transformed, the further evolution is only driven 
by differences in grain boundary energies and becomes much 
slower. As one can see, most grain boundaries are already 
straight (especially high angle grain boundaries between dark 
and light regions). A typical change in structure is shown 
in the two plots at the bottom: snapshots at two successive 
times of the same region (only a small part of the system is 
shown) are displayed. Between the two snapshots, the length 
of the high angle grain boundary (between the dark and light 
region) has significantly decreased, which has led to a motion 
from right to left of the triple junction. 



C. Polycrystals 



We now turn to a more illustrative part, where we show 
results of numerical simulations for the case of isothermal 
solidification (that is, we set u to a constant that does 
not evolve with time) . The aim is to show that our model 
is capable of reproducing the evolution of polycrystals, at 
least qualitatively. Hence, we have simulated the solidifi- 
cation of an undercooled melt with a few seeds of different 
crystalline orientations. The result is shown in Fig. [TO] 
First, the individual grains grow. Since we do not solve 
the heat equation ([5]). no diffusive instability can appear 
and create dendrites, and the shape of the grains remains 
convex. Once the grains impinge, grain boundaries ap- 
pear. They evolve toward straight segments, the lengths 
of which change through the motion of trijunctions while 
the orientation in each grain remains unchanged. Other 
situations such as directional solidification in a fixed tem- 
perature gradient were also simulated, and also led to 
qualitatively correct behaviour of the model. 



IV. CONCLUSION 

We have presented a phase-field model for the solidifi- 
cation and coarsening of polycrystals that is formulated 
in terms of two continuous helds, a phase field that in- 
dicates the local state of matter and an orientation held 
that gives the local direction of the crystallographic axes. 
Contrary to previous orientation-field models found in 
the literature, the free-energy functional of our model 
does not contain a term proportional to the modulus of 
the orientation gradient, but only a standard square gra- 
dient term. Stable localized grain boundary solutions 
are instead generated by a singular coupling function be- 
tween the orientation and the phase field. 

As a result, the mathematical structure of the model 
and the evolution equations are closer to the standard 
phase-field models for solidification, in that the evolution 
equation for the phase field is a regular reaction-diffusion 
equation instead of a singular diffusion equation that has 
to be regularized in a suitable way. The properties of the 
grain boundary solutions can be investigated by standard 
mathematical methods, in particular by the use of the 
well-known mechanical analog that exploits the relation 
between interface solutions and the motion of a particle 
in a potential landscape. This method has allowed us 
to perform a thorough analysis of the grain boundary 
properties, and to choose a suitable mathematical form 
of the singular coupling function. 

One important result of this analysis is that the grain 
boundary solutions of our model have a similar struc- 
ture as the solid-liquid interfaces in standard phase-field 
models: a central region with strong gradients of both 
fields is surrounded by exponential tails, and the inter- 
action between two grain boundaries becomes exponen- 
tially small with their distance. Two important differ- 
ences with previous orientation-field models are a conse- 
quence of these facts. First, the artificial grain rotation 
that is mediated by long-range interactions between grain 
boundaries is not present in our model. Second, it is pos- 
sible to incorporate interfacial anisotropy in our model 
without breaking the variational structure of the model. 
The anisotropy generates unphysical "torque" terms that 
drives the orientation fields towards the minima of the 
surface free energy, but these effects are negligibly small 
for typical simulation parameters that can be achieved in 
practice. 

We have presented some preliminary simulations per- 
formed in two dimensions, which show that a qualitative 
description of multi-grain growth and of the evolution of 
polycrystalline solids can be obtained using this model. 
It should be stressed that the evolution equation for the 
orientation obtained from the variational formalism, Eq. 

is likely to be incorrect, as discussed in more detail in 
Ref. [HI . Therefore, it is unlikely that this model can be 
used for quantitative simulations of grain growth. Nev- 
ertheless, in our opinion this model is interesting because 
of its intrinsic mathematical structure and because of its 
ability to reproduce a host of complex pattern-formation 
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phenomena with the help of a simple set of equations. 

We have concentrated here on the discussion of the 
mathematical structure of the model, but we have neither 
tested its efficiency nor its quantitative accuracy in the 
case of known model problems, such as the growth of an 
isolated dendrite [13] • These are promising subjects for 
further studies. 



Appendix A: Choice of model parameters: further 
details 

The choice of the coupling function for values of the 
phase field close to one (the solid) has already been dis- 
cussed in the main text. Here, we briefly discuss the form 
of the coupling function g close to zero and the double- 
well potential used in our simulations. First, g must be 
equal to zero in the liquid (there should be no energy 
penalty associated with orientation changes in the liquid 
phase). Furthermore, we must have g'(0) — in order 
to keep a minimum of the free energy at the liquid state, 
irrespective of the configuration of the orientation field, 
which has no meaning inside the liquid. Finally, we also 
require that the variation of the orientation should re- 
main confined in the grain boundary. A similar problem 
was solved in a phase-field model of fracture, where the 
strain must be localized inside the crack for a broken 
elastic material. According to the results of Ref. [36j j . 
this imposes that g should behave as (j>@ for (f> — > 0, with 
/3 > 2. Choosing (3 = 2 would lead to a singular be- 
haviour of 9 in the middle of the grain boundary (and to 
strong pinning of the grain boundary on the grid points) . 
Therefore, we choose (3 = 3. 

The conditions on the double- well potential are the fol- 
lowing: it should have two minima at <j) — and <p = 1 
for all temperatures, with a free energy difference pro- 
portional to u between the two wells. Furthermore, its 
second and third derivatives in the solid should be con- 



stant; this is guaranteed if its leading order expansion in 
4> = 1 is of the form (1 — <fi) 2 — 2(1 — <fr) 3 ; with this choice, 
Eq.[23]is satisfied. In addition, we want to prescribe both 
the second derivative of the potential Cn q and its value 
Vo in cf> = 0. The first is set to the same value as the 
second derivative in the solid well, in order to guarantee 
a symmetric interface profile of the phase field, and the 
second is a function of temperature, whereas the depth of 
the solid well is independent of temperature. Assuming 
that V has a polynomial form, the number of equation it 
has to satisfy is 7: 

1. The free energy of the liquid phase is Vq ■ V(0, T) = 
Vo(T) 

2. The liquid phase is a minimum : V^'(0) = 

3. The curvature of the potential at <p — is fixed : 
V"(0) = 2C liq 

4. The free energy of the solid state is 0:1^(1) = 

5. (f> — 1 is a minimum : V'(l) = 

6. The expansion of V in <j> = 1 has to be (1 — <fr) 2 — 
2(1 - <j)f : V"(l) = 2 and = 12 

The first 3 equations impose the value of the constant and 
of the coefficients and 4> 2 in the polynomial expansion of 
V (</)), while the last four equations (the last three points 
of the above enumeration) can be easily translated into a 
set of four independent linear equations on the coefficient 
of V. An evident solution is then to consider that V is 
a 6 th order polynomial and to solve the linear system 
imposed by the last four equations in order to compute 
the remaining coefficients of V. This approach leads to 
the following expression of V: 

V{x) = V + Cnqx 2 + cx 3 + dx 4 + ex 5 + fx 6 . (Al) 

The remaining 4 coefficients have then to satisfy a set 
of 4 independent linear equations that can be inverted 
giving the following polynomial form for V : 



V(x) = U Uq + C\ iq x 2 + cx 3 + dx A + ex 5 + fx 6 . 
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Choosing to have Cn q = 1, we have: 

V{cj>) = V + 4> 2 



(F(1)=0) 
(F'(l) = 0) 
(V"(l) = 2) 

(y( 3 )(i) = 12) 



(A2) 



(2 + 20V )<I> 3 + (1 + 451/ o )0 4 - 36V O 5 + 10V </> 6 



(A3) 



The choice of Vo = —Am (which yields the correct free finally: 
energy difference between the two bulk phases) gives Eq. 

(flUl) . With this choice of V, a reasonable choice for g is 7<j) 3 — 6cA 4 

9 = tt^t (A4) 
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Appendix B: Discretization issues 

Here, we illustrate briefly how a discretization scheme 
that does not derive from an energy can lead to spuri- 
ous motion of the interface while a discretization scheme 
that derives from a free energy cannot lead to any such 
effect. The free energy we consider (for simplicity, in one 
dimension) writes 



T 



dxe(6) 



(V0) 2 



(Bl) 



When simulating the behaviour of a partial differential 
equation (PDE) that derives from this free energy, one 
can either first derive the PDE and then discretize it or 
first discretize the free energy and then derive the evolu- 
tion equation for the discretized field. In the following, 
we consider both options and show how they differ in the 
case of the first term in the above functional. The evo- 
lution equation that derives from this term for the phase 
field 4> is 

d t ct> = d x {e{6)d x <j>) (B2) 
= e'(6)d x 9d x ct> + s(6)dl<j>, (B3) 

where we have used the chain rule to obtain the second 
expression. The latter equation is simply discretized as 
follows: 



d t (j> n = e'(6 n ) 
+ e(8 n )- 



4(Ax) 2 

t>n+l - ^4>n + <j>n-l 



(B4) 



(Ax) 2 

The other approach is to first discretize the free energy: 

2 



E 



2+1/2) / 



Ax 



(B5) 



where 8 n +i/2 is either the value of on a staggered grid 
or the mean value of 9 between the grid points. This, 
using the relaxation law d<fi n /dt = — dTd/d(j) n , leads to 
an evolution equation for <f) n that writes (taking into ac- 
count all the terms of the sum): 



'M+l 



(Ax) 



+ e(0„_ 1 /a)- 



(Ax) 2 

(B6) 

This equation differs from Eq. (|B4j) by little. Indeed, it 
is easy to see that it corresponds to the discretization of 
Eq. (|B2[) instead of Eq. (|B4|) . Using a Taylor expansion 
centered at the grid point n one can show that both dis- 
cretization schemes are equivalent up to the second order, 
but differ by terms of higher order. This just reflects the 



fact that the chain rule is valid for continuous fields, but 
not for discrete fields. 

These seemingly small differences have important con- 
sequences. Due to its construction, Eq. (|B6[) implies that 
Td is decreasing with time: 



= > x 



E 



d<j) n 
2 



dt 



< 0. 



(B7) 
(B8) 



Therefore, the free energy is guaranteed to decrease to- 
ward a minimum. Furthermore, the £ 2 -norm of the rate 
of change of <f> is proportional to the rate of change of 
the discretized energy. This means that if the interface is 
moving, then its squared velocity is proportional to the 
rate of change of the free energy. In other words, any 
motion of the interface (whatever are the discretization 
errors) comes with a decrease in the free energy. This dif- 
fers qualitatively from the expression of Eq. (|B4|) which 
a priori has no reason to derive exactly from any energy 
and where the interface can have a motion proportional 
to discretization errors, as was indeed found in our early 
investigations of this model (35j | . 



Appendix C: Simulations at constant volume of solid 

The discretized equation of motion for the phase field 
has the structure 



d(j) n 
dt 



A. 



B n lL Ti 



(Cl) 



where a one-dimensional example is considered for sim- 
plicity, A n and B n are coefficients that depend on the 
local configuration of phase and orientation fields, and 
U n is the discretized temperature field. For a constant 
but time-dependent temperature, u n = u(t), we have 



dt 



4> n = E An + u W E Bn - 



(C2) 



The left hand side can be interpreted as the change in 
the total volume of solid. Clearly, J? </>„ can be kept 
constant by the following procedure: (i) evaluate the 
coefficients A n and B n , (ii) calculate the temperature 
u = — J2 n 4/(En Bn) which makes the left hand side 
vanish, and (iii) timestep the equations with this value 
of the temperature. After a short time, the temperature 
converges to a constant value, which corresponds to the 
equilibrium temperature for the chosen total volume of 
solid. 
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